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ABSTRACT 

We have studied the effects of gas density inhomogeneities on the escape of ionising 
Lyman continuum (Lyc) photons from Milky Way-type galaxies via 3D numerical 
simulations using the Monte Carlo radiative transfer code CRASH (Ciardi et al. 2001). 
To this aim a comparison between a smooth Gaussian distribution (GDD) and an 
inhomogeneous, fractal one (FDD) has been made with realistic assumptions for the 
ionising stellar sources based on available data in the solar neighborhood. In both 
cases the escape fraction /e^c increases with ionisation rate J\f^ (although for the FDD 
with a flatter slope) and they become equal at — 2 x 10'''" s^^ where fesc — 0.11. 
FDD allows escape fractions of the same order also at lower N^, when Lyc photon 
escape is sharply suppressed by GDD. Values of the escape fraction as high as 0.6 
can be reached (GDD) for jif^ « 9 x 10^" s~^, corresponding to a star formation rate 
(SFR) of roughly 2 Mq yr~^; at this ionising luminosity the FDD is less transparent 
{fesc ~ 0.28). If high redshift galaxies have gas column densities similar to local ones, 
are characterized by such high SFRs and by a predominantly smooth [i.e. turbulence 
free) interstellar medium, our results suggest that they should considerably contribute 
to - and possibly dominate - the cosmic UV background. 
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1 INTRODUCTION 

Recent studies suggest that the ionising radiation escaping 
from galaxies could give a substantial (possibly dominant) 
contribution to the ultraviolet background radiation (UVB) 
during a large cosmological time span (Giallongo, Fontana 
& Madau 1997; Giroux & ShuU 1997; Bianchi, Cristiani & 
Kim 2001). Also, at large redshifts, radiation from the first 
stellar objects has very likely driven the process of cosmic 
reionisation (Gnedin & Ostriker 1998; Ciardi et al. 2000; 
Miralda-Escude, Haehnelt & Rees 2000; Gnedin 2000; Ben- 
son et al. 2000; Ciardi et al. 2001). In spite of this exten- 
sive body of work on both aspects, their predictive power is 
jeopardized by the persisting (theoretical and experimental) 
ignorance on the value of /esc, the fraction of hydrogen- 
ionising photons that escapes from the parent galaxy into 
the intergalactic medium (IGM). Obviously, this quantity 
enters the modeling of the UVB and reionisation process 
and the results depend quite sensibly on the assumptions 
made about this poorly constrained parameter. 
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A promising way to make theoretical progresses on this 
issue is to improve the degree of realism of the modeling and 
the treatment of physical processes to make predictions that 
can be directly compared with available local and interme- 
diate redshift data. 

Dove & Shull (1994b) initially tackled the problem of 
determining fesc by assuming smoothly varying H I distri- 
butions in the Galactic disk. They concluded that about 
14% of Lyman continuum (Lyc) photons are able to escape 
the disk. Dove, Shull & Ferrara (2000) later improved the 
calculation by solving the time-dependent radiation transfer 
problem of stellar radiation through evolving superbubbles. 
Their main result is that the shells of the expanding super- 
bubbles quickly trap ionising photons, so that most of the 
radiation escapes shortly after the formation of the super- 
bubble. This results in a value of fesc roughly a factor of 
2 lower than obtained by Dove & Shull (1994b). Additional 
theoretical works (Ricotti & Shull 2000; Wood & Loeb 2000) 
have extended the analysis to include high-redshift galaxies, 
for which the escape of Lyc photons can be modified by ge- 
ometrical effects {i.e. disk vs. spheroidal systems) and by 
the higher mean galactic interstellar medium (ISM) density. 
Generally speaking, the escape fraction decreases as the ob- 
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Figure 1. Slice of tlie Gaussian (left panel) and fractal (right panel) density distribution described in the text, cut through a plane 
perpendicular to the Galactic disk. The linear size is 1 kpc. Images are displayed with the same logarithmic scale. 



ject virialization redshift or mass become larger; for a typi- 
cal 2-cr fluctuation in a CDM model at redshift ^ 10 these 
studies derive even lower values (/esc ^ 1%). 

Observationally, a wide range of values for the escape 
fraction has been deduced, but it appears that some data 
favor larger values than expected from theory. For example, 
most of the detections of starburst galaxies with the Hopkins 
Ultraviolet Telescope and the Far Ultraviolet Spectroscopic 
Explorer (Leitherer et al. 1995; Hurwitz, Jelinsky & Dixon 
1997; Heckman et al. 2001) are consistent with fesc < 10%, 
although objects have been observed with /esc < 57% (Hur- 
witz, Jelinsky & Dixon 1997). However, absorption from un- 
detected interstellar components could allow the true escape 
fractions to exceed these upper limits. Bland-Hawthorn & 
Maloney (1999; see Erratum 2001) used optical line emis- 
sion data for the Magellanic Stream to derive /esc ~ 45%. 
Even more puzzling are the recent results (Steidel, Pettini 
& Adelberger 2001; Haenhelt et al. 2001) who detected flux 
beyond the Lyman limit (with significant residual flux at 
A < 912A) in a composite spectrum of 29 LBGs ai z — 3.4. 
Also, this implies that at these early epochs galaxies were 
much more transparent to ionising radiation than at present 
time, contrary to the trend found by the theoretical works 
mentioned before. 

Can the low values of fesc predicted by theory and 
the high values suggested by recent observations be brought 
into agreement? In this paper we try to ascertain if the di- 
chotomy can be removed by considering the effects of inho- 
mogeneities in the ISM structure and a more realistic dis- 
tribution of massive stars, the primary sources of ionising 
photons. As already pointed out before, essentially all stud- 
ies (with the partial exception of Wood & Loeb 2000) to 
date have derived the escape fraction assuming smoothly 
varying gas density distributions. However, this position is 
clearly untenable in the light of the large number of observa- 
tions showing that the ISM in the Galaxy and nearby ones 
has a hierarchical, very likely fractal, structure (Elmegreen 



& Falgarone 1996 and references therein). This view is also 
supported by theoretical works (see for example Norman 
& Ferrara 1996), which predict a fractal behavior over a 
large dynamical range. This type of organization naturally 
arises from turbulent pressure, which plays a crucial role 
in the various phases of the ISM (see for example Kulka- 
rni & Fich [1985] for the cold neutral HI component and 
Reynolds [1985] for the warm ionised medium): Norman & 
Ferrara (1996) showed that the ISM turbulent pressure is 
roughly 30 times higher than the thermal one. Turbulence 
is pumped into the ISM primarily by multi-SN explosions; 
however, reaching the high energy density of present day 
galaxies requires a considerable fraction of the Hubble time 
at high redshift. Hence the ISM of these primordial systems, 
differently from their local counterparts, is likely to be more 
quiescent and smooth. 

For these reasons, it seems worthwhile to examine a 
model for the escape fraction in which a fractal, turbulence 
dominated ISM is taken into account. This is what we ex- 
plore in the rest of the paper. 



2 MODEL 

In this section we describe the model adopted for the Milky 
Way (MW) gas density and the stellar distribution, as well 
as the source emission properties. We use the MW as a ref- 
erence template, but our results should hold with good ap- 
proximation for the entire class of MW-type galaxies, as long 
as the properties of their ISM are similar. 

To study the escape of (ionising) Lyman continuum 
(Lyc) photons from the MW, we have chosen a representa- 
tive cubic volume of 1 kpc^ located above the disk midplane 
and centered at the Solar position. This allows us to exploit 
the high quality data available for the spatial distribution of 
the ionising OB stars {e.g. Garmany, Conti & Chiosi 1982). 
Also, the same argument applies to the exquisite detail with 
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which the gas density distribution has been derived (Dickey 
& Lockman 1990; Reynolds & Haffner 2001 and references 
therein). 

Comparison between the soft X-ray background and ab- 
sorption line experiments toward stars near the Sun, have 
suggested that we live in a rarefied, hot ISM cavity, usually 
known as the Local Bubble, whose average radial extent is 
100 pc (Sanders et al. 1977; Snowden et al. 1990; Sfeir et 
al. 1999). As the volume of the Local Bubble is only 0.2% of 
the simulation volume, we neglect for simplicity the effects 
of such structure on our results. 

2.1 ISM Density Distribution 

In order to compare the effects of inhomogeneities on the 
escape of Lyc photons with previous results, we have first 
considered the commonly adopted Gaussian HI density dis- 
tribution (GDD): 

nHi{z) = noexp{-z^ /2H'^), (1) 

where z is the height above the midplane of the disk. 
Such distribution is a solution of hydrostatic equation for 
the gas in the gravitational field of a disk galaxy, includ- 
ing a dark matter component. The Gaussian distribution 
that most closely resembles the three-component Dickey- 
Lockman (1990) vertical distribution has parameters no = 
0.367 cm"^ and H = 0.184 kpc (Dove, ShuU & Ferrara 
2000). This gives a total neutral hydrogen column density 
perpendicular to the disk of Nhi = 5.2 x lO'^" cm~^. 

Observations of molecular clouds have shown that with 
increasing spatial resolution the gas distribution breaks up 
into substructures of yet smaller scales. The self-similarity of 
the relations linking some of the molecular cloud properties 
(such as the mass-size relation and the distribution functions 
of size and mass) suggests that the density distribution very 
likely has a fractal structure induced by turbulence (Falgar- 
one, Phillips & Walker 1991; Elmegreen & Falgarone 1996; 
Elmegreen 1999). As a more realistic description, we have 
thus alternatively adopted a fractal model. 

A fractal ISM distribution is obtained from hierar- 
chically clustered points, along the prescriptions given by 
Elmegreen (1997). In brief, starting from a point with coor- 
dinates (io,io) = (0.5, 0.5), fractals on the x ~ y plane were 
made with K — 6 hierarchical levels by choosing continu- 
ous coordinates such that for random numbers Ri in the 
interval [0,1], ii = io -I- 2{Ri - 0.5)/L^ at the first level, 
i2 ~ ii + 2{R2 — 0.5)/L^ at the second level, and so on up 
to level A"; L = 2 is a geometric factor for subdivision of 
one level into the next. The same procedure was used for 
j. At each level, N = 5 new positions were substituted for 
each position in the previous level, leaving A^^ positions af- 
ter K levels. At each height z, a = 64^ square grid 
was placed around all the points, and the number of points 
inside each grid cell was counted. We thus recovered the 
gas density in each grid cell constraining its average value 
to be the same as for the Gaussian distribution of eq. |l| 
The resulting gas distribution has thus a fractal dimension 
D = logA/logL = 2.32, which is consistent with the experi- 
mentally derived value (Elmegreen & Falgarone 1996). This 
procedure is repeated at each of the 64 heights that have dis- 
crete values in the range < z < 1 kpc; this range represents 
a good compromise between the need to encompass most of 
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Table 1: Simulation Parameters 



Run 




(2) 


E* (^) 


A 


Gaussian 


Gaussian 


24 


B 


Fractal 


Gaussian 


24 


C 


Gaussian 


Gaussian 


48 


D 


Fractal 


Gaussian 


48 


E 


Fractal 


Schmidt Law 


24 



'^'Density Distribution 

(^'stellar Distribution 

(^'stellar Surface Density [kpc~^] 



the vertical extent of the HI disk and to achieve a sufficient 
spatial sampling. The clumping factor of this distribution, 
C = (n^)/ (n)^, calculated as a function of distance from the 
midplane, ranges between 4 and 8. 

In principle, it would have been possible to use the HI 
distribution directly obtained from 21-cm line observations. 
However, resolving power limitations make such maps of lit- 
tle use to study the effects of inhomogeneities. For example, 
the spatial sampling of the Dwingeloo survey (Burton & 
Hartmann 1994) is 0.5° (corresponding to ~ 0.03 kpc) but 
the spectral sampling is 1.03 km s^^, or 0.1 kpc. This is 
about a factor of 6 lower than the resolution required by 
our simulation. Surveys with higher spatial resolution (1 ar- 
cmin) do exist, but they do not cover a sufficiently large 
latitude range (the best is probably the Canadian Galactic 
Plane Survey which covers —3.5° to +5.5°, equivalent to half 
the simulation box) at roughly the same spectral resolution. 
Thus, we have restricted our analysis to the previous two 
model cases. A summary of the adopted density distribu- 
tions in the various simulation runs is shown in Table 1. For 
illustration, in Fig. ^ we show a slice of the GDD (left panel) 
and the FDD (right panel) cut through a plane perpendic- 
ular to the Galactic disk. 

Measurements of faint optical interstellar emission lines 
and pulsar dispersion measures, show that in the MW ap- 
proximately 1/3 of the HI mass is contained in the Warm 
Ionized Medium (Reynolds 1991a; Reynolds 1993). To ac- 
count for the presence of such component, we have multi- 
plied the GDD and FDD described above by a factor of 1.3. 
We assume that initially the gas has a temperature of 100 K, 
as derived for the Cold Neutral Medium. 



2.2 Stellar Distribution 

Garmany, Conti & Chiosi (1982) compiled a catalog of O- 
type stars in the Solar neighborhood. The catalog is sup- 
posed to be complete to a distance of 2.5 kpc from the Sun. 
Within this radius, the surface density of such component 
is about = 24 stars kpc~^. Therefore, we have repro- 
duced the local stellar population by randomly extracting 
positions and photon ionisation rates for 24 stars in our 1 
kpc'' computational volume as follows. 

We have adopted for the stars a vertical Gaussian dis- 
tribution with a scaleheight = 63 pc, as recently derived 
by Mafz-Apellaniz (2001) from a sample of 0-B5 stars from 
the Hipparcos catalog. By randomly sampling this distribu- 
tion, we have derived the z-coordinate of each star; the x and 
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y coordinates have been selected randomly. All stars above 
the galactic plane have been included in the simulations. 

As the star formation process should preferentially take 
place in the densest regions of the ISM, we have also pro- 
duced stellar distributions assuming that the probability of 
finding a star in a certain region is proportional to a power of 
the gas density, P oc p^'^, i.e. a Schmidt-type law (Kennicut 
1998). 

Finally, we have also produced stellar distributions with 
a surface density two times larger than the local one, E* = 
48 kpc"'^, to study the dependence of fesc on a larger range 
of ionisation rates. A summary of the adopted stellar distri- 
butions is shown in Table 1. 



2.3 Source Emission Properties 

A spectral type and luminosity class has been randomly as- 
signed to each star, using the number frequencies of objects 
of a given type in the catalog of Garmany, Conti & Chiosi 
(1982). Within 2.5 kpc from the Sun, about 50% of the stars 
are of O9.0 - 09.5 types, half of which of luminosity class 
V. We have then used the stellar models of Schaerer & de 
Koter (1997) to derive the rate of ionising photons from the 
given spectral type and luminosity class. 

Following this recipe, we have produced 20 different 
stellar distribution realizations for each case in Table 1, to 
get some semblance of the scatter among different realiza- 
tions. The scatter is produced by the random sampling of 
the number, position and ionising photon rate of individual 
stars in the simulation volume. 

For simplicity, we have adopted a single black-body 
spectrum for each star. We have assumed a temperature 
of 40000 K, a mean value for the effective temperatures of 
0-type stars (T^// = 30000-50000 K; Schaerer & de Koter, 
1997). The pattern of ionisation does not heavily depend on 
the detailed shape of the spectrum, as long as most of the 
ionising photons are emitted close to the HI ionisation limit. 

Early B-type stars also produce ionising radiation. As- 
suming a Salpeter Initial Mass Function (IMF) and the ionis- 
ing photon rate of Schaerer & de Koter (1997), we have esti- 
mated that BO-BO. 5 stars contribute on average to less than 
10% of the total number of ionising photons. We have run a 
few simulations to check the influence of B stars. The differ- 
ence between models with and without B stars are masked 
by the random scatter in the results. Therefore, we did not 
include B-type stars in the stellar distributions. 



3 RADIATIVE TRANSFER SIMULATIONS 

To study the propagation of ionising radiation produced by 
the stars through the given ISM density distribution we use 
the Monte Carlo (MC) radiative transfer code CRASH {Cos- 
mological RAdiative transfer Scheme for Hydrodynamics), 
described in Ciardi et al. (2001). The code has been origi- 
nally developed for cosmological applications, but below we 
show that it can correctly handle also the rather different 
conditions prevailing in the interstellar medium of galaxies. 
For clarity, we briefly summarize the main features of the 
numerical scheme relevant to the present study. 

In the application of a MC scheme to radiative transfer 



problems, the radiation intensity is discretized into a rep- 
resentative number of monochromatic photon packets. The 
processes involved (e.g. packet emission and absorption) are 
then treated statistically by randomly sampling the appro- 
priate distribution function. The 1 kpc^ simulation volume 
has been discretized in 64^ grid cells; we consider absorption 
by HI and dust; we neglect the contribution of He because of 
the paucity of the He-ionising photons produced. We have 
estimated the dust optical depth in the optical using the 
canonical relation between the HI column density, Nhi, and 
the color excess E(B-V) (Bohlin, Savage & Drake 1978). This 
has been converted into an optical depth at the ionisation 
limit by adopting the far-UV parametrization of the mean 
Galactic extinction law (Fitzpatrick &Massa 1988). We ob- 
tain: 

rd(912A) = 1.9 X 10"^^(iVH/cm"^), (2) 

where we have excluded the contribution of dust scattering 
to the opacity, by multiplying for {1 — uj)\ uj = 0.4 is the 
average albedo of dust in the UV close to the ionisation 
limit (Witt et al. 1993). Given the above equation, as the 
final mean ionisation fraction for run A (run B) is (i) = 0.69 
(0.40) (see Section 4), for the frequencies of interest here, we 
find that Td/Tui < 2.46 x 10"^ (1.28 x 10"^). Thus, dust 
contribution to absorption is negligible. 

For the range of densities considered here, a number of 
photon packets equal to Afp = 5 x 10^ has to be emitted in 
order to reach numerical convergence. 

Differently from Ciardi et al. (2001), here we deal with 
more than one ionising source. We have found that, if the 
number of emitted packets is sufficiently high, a more re- 
liable solution of the discretized time-dependent ionisation 
equation is obtained if the time step At used in their eq. 10 
is not treated statistically. Thus, here At is the time elapsed 
since a photon packet has gone through the cell for which 
the equation is being solved. Only photons escaping from the 
top side of the cubic region contribute to the escape fraction, 
while those escaping from side faces or the bottom one are 
subject to a reflecting boundary condition, thus simulating 
photons coming from adjacent regions. 

The same numerical tests of the code described in Ciardi 
et al. (2001) have been performed in an interstellar (rather 
than intergalactic) environment. We do not report them 
here; instead, in Fig. |^ we show the evolution of ionised 
regions (in black) produced in the GDD described in the 
previous Section, by sources with different ionisation rate, 
located on the disk midplane. The slices show a plane per- 
pendicular to the Galactic disk through each source location. 
Different columns refer to different source ionisation rates 
(TV't = 1.8, 16, 34 X 10"** s"^ respectively, from left to right), 
while from the top to the bottom we show the ISM ionisation 
structure at times t — 0.05, 1, 1.5, 5 Myr after the source has 
been turned on. The solid lines correspond to the stationary 
analytical solution given in Dove & ShuU (1994b). The ion- 
isation rates were chosen in order to reproduce the various 
possible shapes of the ionised regions: spherical, elongated 
and funnel-like. 

To prevent misunderstandings, it is useful here to give 
our working definition of the escape fraction. We define the 
global escape fraction, fssc,g, as: 
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Figure 2. Evolution of ionised regions (black) produced in the GDD described in Section 2, by sources with different ionisation rate, 
located on the midplane of the disk. The linear box size is 1 kpc. Columns, from the left to the right, refer to a source ionisation rate of 
J\f~^ = (1.8, 16, 34) X 10''® s~^ respectively; from the top to the bottom the ionisation structure at a time t = (0,05, 1, 1,5, 5) Myr after 
source turn on is shown. The solid lines correspond to the stationary analytical solution (Dove & Shull 1994b). 
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fesc,g (t) 



(3) 



the number of 



where TV"^ is the escaping photon rate 
photons reaching 2 > 1 kpc per unit time) and f\f^ is the 
total photon rate production at (simulation) time t. The 
instantaneous escape fraction, feiac,i, is instead given by: 



fesc,i (^) 



M^{t)' 



(4) 



In order to derive fesc, we assume that the stars have 
a constant ionisation rate and we run the simulation un- 
til convergence (i.e. the escape fraction and the ionisation 
structure do not change) is reached. As we find that, after 
a time of ~ 5 x 10® yr the value of f^sc remains roughly 
constant, we stop the simulations after ~ 10^ yr. This is 
expected as the recombination time at the midplane is 0.3 
Myr. 

Following the method described above, we have run 20 
simulations for each of the cases indicated in Table 1. 



4 RESULTS 

In Fig. ^ we show the evolution of the global escape frac- 
tion, /esc.g, as a function of the total ionisation rate, TV^, 
for runs A (filled triangles in Fig. |^) and B (filled circles 
in Fig. ^b). In both cases /esc,g increases with increasing 
J\f~j, although for the FDD with a flatter slope. The reason 
for this different behavior is that photons in a GDD can es- 
cape only if an ionised channel can be produced by the stars 
themselves. Consequently, a low ionisation rate results in a 
low escape fraction. In a FDD instead, photons can travel 
along clear sigthlines through low density ionised channels 



and escape more easily. Thus, for 7V-y 



always higher for a FDD, while for larger ionisation rates a 
GDD becomes more transparent. In both cases, the scatter 
in /esc,g is mainly due to differences in the star position and 
ionisation rate. 

To study the effect of a larger total ionisation rate, we 
have increased the stellar surface density, as explained in 
Section 2, obtaining a maximum value M-y ~ 5 x lO^'^ s^^. 
These additional runs are represented in Fig. ^ as open tri- 
angles (run C) and circles (run D). The scatter is reduced 
with respect to runs A and B as, increasing the number of 
stars, their position affects less sensibly the value of fesc, 
which instead depends more strongly on the total ionisation 
rate. In order to check if the curve for the GDD would reach 
a plateau or rather keep growing (as the analytical deriva- 
tion by Dove & ShuU 1994a suggests), we have increased the 
total ionisation rate up to a value of ~ 8 x 10^" s~^ (cor- 
responding to a stellar surface density of ~ 120 kpc~^). 
As we expect the scatter to be reduced, in order to limit 
the computational time, we have only run 6 different real- 
izations. 

The evolution of /esc, 9 is linear over the considered 
range of jif-^ and it can be fitted by the following function: 



(5) 



where for a GDD (FDD) a = 0.089 ± 0.003 (0.023 ± 0.003) 
and (3 = -0.07 ± 0.01 (0.07 ± 0.01). In both cases the fit is 
represented as a dotted line in Fig. M. We find that the GDD 
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Figure 3. (a) Evolution of the global escape fraction, fesc,g, 
as a function of the total ionisation rate, Af-y, for runs A (filled 
triangles) and C (open triangles). The dotted line is the fit given 
by eq. |^. The dashed line is the analytical solution derived by 
Dove &; ShuU (1994a), assuming a single source located on the 
galactic plane, (b) Same as (a) for runs B (filled circles) and D 
(open circles). 



curve eventually flattens and departs from the flt in eq. |5| 
on the contrary the results for the FDD are well described 
by the same flt also for high values of jV^. 

Dove & ShuU (1994a) derived analytically the value of 
/esc expected in a GDD from a point source of given ionisa- 
tion rate, located on the disk of the Galaxy. The value they 
obtained (represented by a dashed line in Fig. ^) should 
be compared with the asymptotic value of /esc,i. However, 
as for late evolutionary times /esc, 9 approaches /esc,i (see 
below), we can compare their results with our runs A and 
C directly from Fig. ^a. Clearly, the evolution of /esc, 9 dif- 
fers from the analytic calculation and higher values for the 
escape fraction are obtained, especially for large ionisation 



© 2001 RAS, MNRAS 000, 



Lyman Continuum Escape 7 



-0.5 



i -1.5 



o 



-2.5 




10^ 



10^ 106 
Time [yr] 



10^ 



Figure 4. Time evolution of the escape fraction for a Gaussian 
(lower curves) and a fractal (upper curves) density distribution. 
The evolution of /esc,g (solid lines), as well as the one of fesc,i 
(dotted lines) is shown. 



rates. This is due to the fact that in our case the ionisation 
rate is distributed among various sources located at different 
heights above the midplane. 

In Fig. ^ the time evolution of the escape fraction is 
shown, for a Gaussian (lower curves) and a fractal (upper 
curves) density distribution. As a reference, we have chosen a 
run that produces, for both distributions and given the same 
set of stars (with Afj = 2 x 10^° s~^), a comparable value 
of fesc,g- The solid (dotted) lines indicate the evolution of 
fesc,g {fesc,i)- As in a GDD photons can escape only through 
ionised channels produced by the stars themselves, the pho- 
ton escape is retarded and the evolution of fesc is steeper 
compared to the case of a FDD. Initially, fesc,i, which be- 
comes constant after ~ 10^ yr, is higher than fesc,g- This 
is due to the fact that, when the first photons escape (at a 
time t ~ 10^ [lO'^] yr for a FDD [GDD]), a large number of 
photons has already been emitted and this results in a very 
low value of fesc,g- Then fesc,g increases, approaching /esc,i. 

In Figs. I and | illustrative slices extracted from the 
simulation box for the same run of Fig. ^, show the evolution 
of the ionised gas (black regions) produced by the sources 
in the given density field. The top six panels are for a GDD 
and show the ionised regions at a time t = 0.1 Myr from 
the source turn on, at different vertical distances from the 
midplane. The bottom six panels are the same for a FDD. 
As already pointed out above, although in the GDD the 
upper gas layers at this time are still neutral, some ionising 
photons have already escaped the FDD. In Fig. ^ the top 
(bottom) six panels refer to the GDD (FDD) . Here we show 
the time evolution (columns, from left to right, refer to t = 
0.005, 0.5, 5 Myr respectively) of slices at two heights above 



the midplane. Again, the ionised regions are more extended 
in a fractal than in a Gaussian medium. 

Although in this work we focus on the variation of 
fesc with the different density distribution, it is interest- 
ing to check the values we obtain for the column den- 
sity of the neutral gas. For runs A and B, we should be 
able to recover the observed Galactic HI column density 
(Nhi = 5.2 X 10^° cm"^ i.e. about 30% of the total gas 



should be ionised (Sect. 2.1). For a GDD (run A) we obtain a 



final mean HI column density of Nhi = 2.1 x 10 cm , i.e. 
about 40% lower; for a FDD (run B), Nhi = 4.1 x 10^" cm"^ 
is closer to the observed value. In both cases, however, some 
of the simulations have column density closer to the ob- 
served. Despite the larger escape fraction for a FDD, the 
resulting mean HI column density is larger for a fractal 
medium, where the highest density regions are more diffi- 
cult to ionise and recombine faster. Also with an increased 
stellar surface density (runs C and D) the mean Nhi in a 
FDD remains as high as « 2.7 x 10^" cm"^, while in a GDD 
it is ~ 5.9 x 10^^ cm~^, with values as low as ~ 10^^ cm~^. 
The reduced ionisation fraction in the case of an inhomoge- 
neous medium was also noted in the model for the diffused 
ionised medium of Miller & Cox (1993), which took into ac- 
count the absorption of UV radiation by clouds with density 
larger than the surrounding medium in a statistical way. 

In Fig. ^ the vertical distribution of the mean HI (solid 
lines) and HII (dotted lines) number densities is shown, for 
runs A and B with the above Nhi, in the case of a Gaus- 
sian (top panel) and fractal (bottom panel) density field at 
different times after the source turn on. For both distribu- 
tions, the HI profile does not change until t « 10^ yr, when 
the inner regions, where the stars reside, have an HII den- 
sity of nnii ~ 0.1 cm~^. The final HI density profile has 
a similar shape for both distributions, with a slightly lower 
value in the inner regions for a GDD case. The HII number 
density increases gradually in a FDD, while in a GDD, as 
previously mentioned, the ionisation of the upper gas lay- 
ers is postponed until the stars have been able to produce 
an ionised channel for the escaping of photons. At the final 
stages, the HII distribution appears more extended than the 
HI one. This is consistent with observations of the Reynolds 
layer, which seems to have a vertical scale height of « 1 kpc 
(Mezger 1978; Reynolds 1991a, 1991b). However, a proper 
modelling of the diffuse ionised gas should take into account 
the dynamical effects on gas, which would make the ionised 
gas to expand, as well as a larger vertical size for the simu- 
lation box. 

In Fig. ^ we show the evolution of fesc,g as a function 
of the total ionisation rate for run B (filled circles) and E 
(stars), the only difference being the stellar distribution. The 
distribution of the escape fraction looks similar in the two 
cases and the obtained values are comparable; in fact, in 
run E the stars, following the Gaussian distribution of the 
gas, have a larger vertical extension, but, as they are located 
in denser regions, more photons are needed to ionised them 
and the two effects roughly balance. 

Finally, we present in Fig. ^ the spectrum of the emerg- 
ing ionising radiation, for five realizations of runs A and B, 
respectively. In both plots, the uppermost curve is the in- 
trinsic spectrum of stars, i.e. a black body with T=40000 K. 
As already discussed, each simulation is characterised by the 
rate of emitted ionising photons, M-,. Since we are mainly 
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Figure 5. Evolution of ionisation field (black) produced by a given stellar distribution into a Gaussian (top six panels) and a fractal 
(bottom six panels) density field. The slices are taken at a time t = 0.1 Myr from the source turn on, at different vertical distances from 
the midplane (0, 10, 156, 312, 625, 940 pc, left to right). 



interested in the shape of the spectrum, we normalise all 
simulations to the value J\fj — 2 x lO^^s"^. As expected 
from the v'"^ frequency dependence of the HI photoionisa- 
tion cross section, the emerging filtered spectrum is harder 
than the intrinsic one. For run A, the shape of the spec- 



trum depends on the amount of escaping photons (i.e. on 
the escape fraction) in each realization: the lower and upper 
emerging spectra in the left panel of Fig. ^ being those for 
the lower and higher values of 7V7. Different realizations of 
run B yield much more similar emerging spectra, this re- 
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Figure 6. Evolution of ionisation field (black) produced by a given stellar distribution into a Gaussian (top six panels) and a fractal 
(bottom six panels) density field. The time evolution (columns, from left to right, refer to t = 0.005, 0.5, 5 Myr respectively) of slices at 
(first row) and 265 pc (second row) above the midplane is shown. 



fleeting the nearly constant value of the escape fraction. As 
a reference, we plot in the right panel of Fig. ^ a black- 
body with T=45000 K. We remind that here we do not take 
into account the absorption of photons with hi/ > 24.6 eV 
by Hel, which may be the cause for the low Hel ionisation 



observed in the diffuse ionised medium (Reynolds & Tufte 
1995). 
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Figure 9. Spectrum of the emerging ionising radiation, for five realizations of run A (left panel) and run B (rigth panel). In both panels, 
the uppermost curve is the spectrum of the intrinsic ionising radiation (a blackbody with T=40000 K). For ease of presentation, all the 
spectra have been normalised to an intrinsic photon rate J\f^ = 2 X lO^'^s"^. The dotted curve in the right panel is a blackbody with 
T=45000 K. 
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Figure 7. Vertical distribution of the mean HI (solid lines) and 
HII (dotted lines) number densities, in the case of a Gaussian 
(top panel) and fractal (bottom panel) density field. The numbers 
refer to different times from the source turn on: i = 0...6 refers, 
respectively, to i = 0, 10^, 10^, 10*, 10^, 10*' and lO'^ yr. 



5 SUMMARY AND CONCLUSIONS 

We have studied the eflfects of gas density inhomogeneities 
on the escape of ionizing Lyman continuum photons from 
Milky Way-type galaxies via radiative transfer numerical 
simulations. To this aim a comparison between a smoothly 
stratified and an inhomogeneous, fractal distribution have 
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Figure 8. Evolution of the global escape fraction, fesc,g, as a 
function of the total ionisation rate, Af-y, for runs B (filled circles) 
and E (stars). 



been made with realistic assumptions for the ionizing stellar 
sources based on available data in the solar neighborhood. 
The main results obtained can be summarized as follows. 

(1) The global escape fraction, fe3c,g, in the case of a 
Gaussian Density Distribution (GDD) rapidly increases with 
increasing total ionization rate, A/'^., and it flattens for high 
values of A/'^; for a Fractal Density Distribution (FDD) the 
dependence on A''^ is milder. 

(2) For Af-y ^ 2 X 10^° s~^ fesc,g is always higher for a 
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FDD, while for larger ionization rate a GDD becomes more 
transparent. 

(3) For low stellar surface densities, E*, there is a large 
scatter in the results, due to differences in the star positions 
and ionization rates. The scatter is reduced for larger E^, as 
the star position affects less sensibly the value of /esc, which 
instead depends more strongly on J\f^. 

(4) In a GDD the photon escape is retarded and the 
time evolution of the escape fraction is steeper compared to 
the FDD case. 

(5) For a GDD we obtain a final mean HI column den- 
sity of Nhi = 2.1 X 10^" cm~^, lower than the observed one; 
for a FDD, Nhi = 4.1 x 10^" cm~^ closely matches the ex- 
perimental value. In all cases, the HII distribution appears 
to be more extended than the HI one. 

(6) The value of the escape fraction is more sensitive to 
the gas density rather than the stellar distribution. 

In order to draw a more general conclusion from the 
results summarized above, it is useful to translate the to- 
tal ionization rate into a star formation rate for a given 
galaxy. To do so we use the results of the stellar population 
synthesis code Starburst99 (Leitherer et al. 1999). Assum- 
ing a continuous star formation mode, a Salpeter IMF with 
Mdown ~ IMq and Mup = IOOMq and a metallicity equal 
to 1/20 of solar, we find that the expected ionization rate 
inside our simulated kpc'' produced by the corresponding 
stellar population is Af^ = 4.5 x 1O^°(M*/M0yr"^). From 
Fig. ^ we then note that for M* « 2MQyr~^ the escape 
fraction is already equal to roughly 0.6 (0.28) for the GDD 
(FDD) case. 

Star formation rates of a few tens of solar masses are 
commonly derived from observations of LBGs (Pettini et 
al. 2001). This result may suggest a galaxy-dominated UV 
background. In fact, Bianchi, Cristiani & Kim (2001) have 
recently shown that estimates of the local and high-z meta- 
galactic ionizing flux are consistent with a galaxy-dominated 
background if feac ~ 10%. Hence it appears that these type 
of sources, as long as they could be modelled as relatively 
normal disk galaxies, can strongly influence the ionization 
history of the intergalactic medium and, possibly, the galaxy 
formation process. 

In principle, then, photons can escape rather efficiently 
(particularly if their ISM is smooth) from high redshift 
galaxies, provided their vertical gas distribution is not too 
different from the Milky-Way type galaxies investigated 
here. Although galactic disk models based on simple semi- 
analytical prescriptions (Fall & Efstathiou 1980; Mo, Mao 
& White 1998; Weil, Eke & Efstathiou 1998; Ferrara, Pet- 
tini & Shchekinov 2000; Wood & Loeb 2000) predict that 
the disk column density should increase with redshift by a 
factor ~ 10 between the present value and redshift 5, this 
estimate is still subject to several uncertainties. For exam- 
ple it is not clear how both the gas temperature (governing 
the disk scaleheight) and the fraction of baryons that are 
able to cool and settle in the disk (determining the mid- 
plane gas density) evolve with cosmic time. Even if the col- 
umn density is larger, we do expect that the galactic ISM is 
smoother than the present day one, due to the short time 
interval available for the build up of a fully developed tur- 
bulent spectrum and hence a fractal gas distribution. This 
process might require several (say, 10) eddy turnover times, 
te, on the scale at which turbulence is injected by super- 



bubbles, i.e. £ ^ 1 kpc (Norman & Ferrara 1996). If the 
characteristic eddy velocity is iif ~ 1 km s~^, this implies 
that te ~ 10 Gyr is much longer than the Hubble time at 
redshift 3. Hence, a final conclusion from models like the one 
presented here extended to encompass high redshift galaxies 
have to await additional data (both observational and from 
galaxy formation models) on the distribution and dynamical 
state of the gas in primordial galaxies. 
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